home *** CD-ROM | disk | FTP | other *** search
/ Space & Astronomy / Space and Astronomy (October 1993).iso / mac / programs / space / DE118I.ZIP / EPSILN.C < prev    next >
C/C++ Source or Header  |  1993-02-06  |  7KB  |  254 lines
  1. /* Obliquity of the ecliptic at Julian date J
  2.  *
  3.  * IAU Coefficients are from:
  4.  * J. H. Lieske, T. Lederle, W. Fricke, and B. Morando,
  5.  * "Expressions for the Precession Quantities Based upon the IAU
  6.  * (1976) System of Astronomical Constants,"  Astronomy and Astrophysics
  7.  * 58, 1-16 (1977).
  8.  *
  9.  * Before or after 200 years from J2000, the formula used is from:
  10.  * J. Laskar, "Secular terms of classical planetary theories
  11.  * using the results of general theory," Astronomy and Astrophysics
  12.  * 157, 59070 (1986).
  13.  *
  14.  *  See precess.c and page B18 of the Astronomical Almanac.
  15.  * 
  16.  */
  17.  
  18. #ifndef NOINCS
  19. #include "mconf.h"
  20. #include "prec.h"
  21. #include "const.h"
  22. #endif
  23.  
  24. #if LDOUBLE
  25. #if IBMPC
  26. short epsAi[] = {
  27. 0xb229,0x50d6,0x2f6a,0xeda2,0x3ff5, XPD
  28. 0xedd0,0x8d25,0x3ad1,0x9aaa,0xbff4, XPD
  29. 0xc28f,0x28f5,0x8f5c,0xbb42,0xc004, XPD
  30. 0x4dd3,0x1062,0xb958,0xa4ce,0x400f, XPD
  31. };
  32. short epsNi[] = {
  33. 0xfbf4,0xcdfe,0x138b,0xed5f,0x3ff5, XPD
  34. 0xedd0,0x8d25,0x3ad1,0x9aaa,0xbff4, XPD
  35. 0xe148,0x147a,0x47ae,0xbb61,0xc004, XPD
  36. 0xa4a9,0x404e,0x2113,0xa4e6,0x400f, XPD
  37. };
  38. short epsLi[] = {
  39. 0x12eb,0x0788,0xaf45,0x86b0,0x3fdf, XPD
  40. 0xf75f,0xd639,0x60eb,0xc6f1,0x3fe3, XPD
  41. 0x43e0,0x4bdb,0x3c80,0x95a0,0x3fe9, XPD
  42. 0x3646,0xb013,0x4476,0xbf20,0x3fea, XPD
  43. 0x7f66,0x2345,0x9e44,0xa3c9,0xbff0, XPD
  44. 0xad85,0x117e,0xacd9,0xa39f,0xbff6, XPD
  45. 0xde1a,0xc1ac,0xaafb,0xa85c,0xbff7, XPD
  46. 0x8106,0x4395,0x6c8b,0xffe7,0x3fff, XPD
  47. 0xc083,0xa1ca,0xb645,0xfdf3,0xbff8, XPD
  48. 0x9581,0x8b43,0xe76c,0xea0b,0xc007, XPD
  49. 0x4dd3,0x1062,0xb958,0xa4ce,0x400f, XPD
  50. };
  51. #endif
  52. #if MIEEE
  53. long epsAi[] = {
  54. 0x3ff50000,0xeda22f6a,0x50d6b229,
  55. 0xbff40000,0x9aaa3ad1,0x8d25edd0,
  56. 0xc0040000,0xbb428f5c,0x28f5c28f,
  57. 0x400f0000,0xa4ceb958,0x10624dd3,
  58. };
  59. long epsNi[] = {
  60. 0x3ff50000,0xed5f138b,0xcdfefbf4,
  61. 0xbff40000,0x9aaa3ad1,0x8d25edd0,
  62. 0xc0040000,0xbb6147ae,0x147ae148,
  63. 0x400f0000,0xa4e62113,0x404ea4a9,
  64. };
  65. long epsLi[] = {
  66. 0x3fdf0000,0x86b0af45,0x078812eb,
  67. 0x3fe30000,0xc6f160eb,0xd639f75f,
  68. 0x3fe90000,0x95a03c80,0x4bdb43e0,
  69. 0x3fea0000,0xbf204476,0xb0133646,
  70. 0xbff00000,0xa3c99e44,0x23457f66,
  71. 0xbff60000,0xa39facd9,0x117ead85,
  72. 0xbff70000,0xa85caafb,0xc1acde1a,
  73. 0x3fff0000,0xffe76c8b,0x43958106,
  74. 0xbff80000,0xfdf3b645,0xa1cac083,
  75. 0xc0070000,0xea0be76c,0x8b439581,
  76. 0x400f0000,0xa4ceb958,0x10624dd3,
  77. };
  78. #endif
  79. #define epsA ((long double *)epsAi)
  80. #define epsN ((long double *)epsNi)
  81. #define epsL ((long double *)epsLi)
  82. #else /* not LDOUBLE: */
  83. #if DEC
  84. static short epsAi[4*4] = {
  85. 0035755,0121057,0065120,0153262,
  86. 0135432,0125072,0150615,0022756,
  87. 0141473,0041217,0056050,0172703,
  88. 0044244,0147271,0054020,0061116,
  89. };
  90. static short epsNi[4*4] = {
  91. 0035755,0057423,0105715,0177374,
  92. 0135432,0125072,0150615,0022756,
  93. 0141473,0060507,0127024,0075341,
  94. 0044244,0163041,0011500,0047245,
  95. };
  96. static short epsLi[11*4] = {
  97. 0030206,0130257,0042407,0104023,
  98. 0031306,0170540,0165726,0034767,
  99. 0032625,0120074,0100113,0155504,
  100. 0033077,0020104,0073260,0011466,
  101. 0134443,0144636,0042043,0042577,
  102. 0136043,0117654,0154421,0077256,
  103. 0136250,0056252,0175701,0126336,
  104. 0040377,0163554,0105503,0112601,
  105. 0136575,0171666,0042641,0145301,
  106. 0142352,0005747,0066213,0041626,
  107. 0044244,0147271,0054020,0061116,
  108. };
  109. #endif /* double DEC */
  110. #if IBMPC
  111. static short epsAi[4*4] = {
  112. 0x1ad6,0xed4a,0xb445,0x3f5d,
  113. 0xa4be,0x5a31,0x5547,0xbf43,
  114. 0x1eb8,0xeb85,0x6851,0xc047,
  115. 0x0c4a,0x2b02,0x99d7,0x40f4,
  116. };
  117. static short epsNi[4*4] = {
  118. 0xbfdf,0x7179,0xabe2,0x3f5d,
  119. 0xa4be,0x5a31,0x5547,0xbf43,
  120. 0x8f5c,0xf5c2,0x6c28,0xc047,
  121. 0x09d5,0x2268,0x9cc4,0x40f4,
  122. };
  123. static short epsLi[11*4] = {
  124. 0xf102,0xe8a0,0xd615,0x3df0,
  125. 0xc73f,0x1d7a,0xde2c,0x3e38,
  126. 0x7b68,0x9009,0xb407,0x3e92,
  127. 0x0267,0x8ed6,0xe408,0x3ea7,
  128. 0x68b0,0xc884,0x7933,0xbf04,
  129. 0x2fd6,0x9b22,0x73f5,0xbf64,
  130. 0x359c,0x5f78,0x0b95,0xbf75,
  131. 0x72b0,0x9168,0xfced,0x3fff,
  132. 0x3958,0xc8b4,0xbe76,0xbf8f,
  133. 0x6873,0xed91,0x417c,0xc07d,
  134. 0x0c4a,0x2b02,0x99d7,0x40f4,
  135. };
  136. #endif /* double IBMPC */
  137. #if MIEEE
  138. static short epsAi[4*4] = {
  139. 0x3f5d,0xb445,0xed4a,0x1ad6,
  140. 0xbf43,0x5547,0x5a31,0xa4be,
  141. 0xc047,0x6851,0xeb85,0x1eb8,
  142. 0x40f4,0x99d7,0x2b02,0x0c4a,
  143. };
  144. static short epsNi[4*4] = {
  145. 0x3f5d,0xabe2,0x7179,0xbfdf,
  146. 0xbf43,0x5547,0x5a31,0xa4be,
  147. 0xc047,0x6c28,0xf5c2,0x8f5c,
  148. 0x40f4,0x9cc4,0x2268,0x09d5,
  149. };
  150. static short epsLi[11*4] = {
  151. 0x3df0,0xd615,0xe8a0,0xf102,
  152. 0x3e38,0xde2c,0x1d7a,0xc73f,
  153. 0x3e92,0xb407,0x9009,0x7b68,
  154. 0x3ea7,0xe408,0x8ed6,0x0267,
  155. 0xbf04,0x7933,0xc884,0x68b0,
  156. 0xbf64,0x73f5,0x9b22,0x2fd6,
  157. 0xbf75,0x0b95,0x5f78,0x359c,
  158. 0x3fff,0xfced,0x9168,0x72b0,
  159. 0xbf8f,0xbe76,0xc8b4,0x3958,
  160. 0xc07d,0x417c,0xed91,0x6873,
  161. 0x40f4,0x99d7,0x2b02,0x0c4a,
  162. };
  163. #endif /* double MIEEE */
  164. #if UNK
  165. DOUBLE epsA[] = {
  166.  1.81300000000000000000E-3,
  167. -5.90000000000000000000E-4,
  168. -4.68150000000000000000E1,
  169.  8.43814480000000000000E4,
  170. };
  171. DOUBLE epsN[] = {
  172.  1.81100000000000000000E-3,
  173. -5.90000000000000000000E-4,
  174. -4.68450000000000000000E1,
  175.  8.44282584000000000000E4,
  176. };
  177. DOUBLE epsL[] = {
  178.  2.45000000000000000000E-10,
  179.  5.79000000000000000000E-9,
  180.  2.78700000000000000000E-7,
  181.  7.12000000000000000000E-7,
  182. -3.90500000000000000000E-5,
  183. -2.49670000000000000000E-3,
  184. -5.13800000000000000000E-3,
  185.  1.99925000000000000000E0,
  186. -1.55000000000000000000E-2,
  187. -4.68093000000000000000E2,
  188.  8.43814480000000000000E4,
  189. };
  190. #else /* not UNK: */
  191. #define epsA ((double *)epsAi)
  192. #define epsN ((double *)epsNi)
  193. #define epsL ((double *)epsLi)
  194. #endif /* UNK */
  195. #endif /* not LDOUBLE */
  196.  
  197. /* The results of the program are returned in these
  198.  * global variables:
  199.  */
  200. DOUBLE jdeps = 0.0; /* Date for which obliquity was last computed */
  201. DOUBLE eps = 0.0; /* The computed obliquity in radians */
  202. DOUBLE coseps = 0.0; /* Cosine of the obliquity */
  203. DOUBLE sineps = 0.0; /* Sine of the obliquity */
  204. extern DOUBLE eps, coseps, sineps, STR, J1900, J2000, Jcentury, Jmillenium;
  205. DOUBLE SIN(), COS(), FABS();
  206.  
  207. epsiln(J)
  208. DOUBLE J; /* Julian date input */
  209. {
  210. DOUBLE T;
  211.  
  212. if( J == jdeps )
  213.     goto done;
  214.  
  215. #if 0
  216. T = (J - J2000)/Jcentury;
  217.  
  218. /* This expansion is from the AA.
  219.  * Note the official 1976 IAU number is 23d 26' 21.448", but
  220.  * the JPL numerical integration found 21.4119".
  221.  */
  222. if( FABS(T) < Two )
  223.     eps = (((epsA[0]*T + epsA[1])*T +epsA[2])*T + epsA[3])*STR;
  224. #endif
  225.  
  226. #if 0
  227. T = (J - J1900)/Jcentury;
  228. if( FABS(T) < Two )
  229.     eps = (((epsN[0]*T + epsN[1])*T + epsN[2])*T + epsN[3])*STR;
  230. /* This expansion is from Laskar, cited above.
  231.  * Bretagnon and Simon say, in Planetary Programs and Tables, that it
  232.  * is accurate to 0.1" over a span of 6000 years. Laskar estimates the
  233.  * precision to be 0.01" after 1000 years and a few seconds of arc
  234.  * after 10000 years.
  235.  */
  236. else
  237. #endif
  238.  
  239. #if 1
  240.     {
  241.     T = (J - J2000)/Jmillenium;
  242.     eps = ((((((((( epsL[0]*T + epsL[1])*T + epsL[2])*T
  243.     + epsL[3])*T + epsL[4])*T + epsL[5])*T + epsL[6])*T
  244.     + epsL[7])*T + epsL[8])*T + epsL[9])*T + epsL[10];
  245.     eps *= STR;
  246.     }
  247. #endif
  248.  
  249. coseps = COS( eps );
  250. sineps = SIN( eps );
  251. jdeps = J;
  252. done: ;
  253. }
  254.